#include <pqxx/pqxx>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <dfsd.h>
#include <vector>
#include <mfhdf.h>
#include <hdf.h>
#include <gdal.h>

#include <sys/types.h>
#include "libpq-fe.h"

#include <iostream>
#include <sstream>
#include <fstream>

#include <ctime>

using namespace std;
using namespace pqxx;

typedef struct
{
    float32 lat;
    float32 lon;

    uint8 land_cover_type1; //  land_cover_type_1
/*
    8-bit unsigned integer,    2400 x 2400
    Number of attributes = 22
        long_name = Land_Cover_Type_1
        units = class number
        valid_range = 0,254
        _FillValue = 255

        water                               = 0     like ocean/sea
        evergreen needleleaf forest         = 1
        evergreen broadleaf forest          = 2
        deciduous needleleaf forest         = 3
        deciduous broadleaf forest          = 4
        mixed forests                       = 5
        closed shrubland                    = 6
        open shrublands                     = 7
        woody savannas                      = 8
        savannas                            = 9
        grasslands                          = 10
        permanent wetlands                  = 11    seem like river/lake
        croplands                           = 12
        urban and built-up                  = 13
        cropland/natural vegetation mosaic  = 14
        snow and ice                        = 15    ? why
        barren or sparsely vegetated        = 16
        unclassified                        = 254   !found
*/
    uint8 land_cover_type2;
    uint8 land_cover_type3;
    uint8 land_cover_type4;
    uint8 land_cover_type5;
    uint8 land_cover_type1_assessment;
    uint8 land_cover_type2_assessment;
    uint8 land_cover_type3_assessment;
    uint8 land_cover_type4_assessment;
    uint8 land_cover_type5_assessment;
    uint8 land_cover_type_qc;
    uint8 land_cover_type1_secondary;
    uint8 land_cover_type1_secondary_percent;
    uint8 lc_property_1;
    uint8 lc_property_2;
    uint8 lc_property_3;


} CombiData;

GByte getHDFMODPixelValueT(char *fMODISName, char *bName, int bNum, VOIDP outRaster);
void    findDimension(char *nameInputFile, int **arrayDimension, char *nameBand);
void    getFactor(char * inputFileName, char * band, char *attribute, float **valueToSearch);
bool    postgresql_query(string);

int main()//int argc)//, char* argv[])
{

    time_t t1 = time(0);    
    int left = 2232;
    int right = 2394;
    int top = 2225;
    int bottom = 2334;

    char        *fMCD12Q1Name = NULL;
    int         *dimDataset, dx, dy, k, j, pos; 
  
    uint8 *land_cover_type1 = NULL, *land_cover_type2 = NULL, *land_cover_type3 = NULL, *land_cover_type4 = NULL, *land_cover_type5 = NULL, *land_cover_type1_assessment = NULL, *land_cover_type2_assessment = NULL, *land_cover_type3_assessment = NULL, *land_cover_type4_assessment = NULL, *land_cover_type5_assessment = NULL, *land_cover_type_qc = NULL, *land_cover_type1_secondary = NULL, *land_cover_type1_secondary_percent = NULL, *lc_property_1 = NULL, *lc_property_2 = NULL, *lc_property_3 = NULL;
    
    GByte           resHdf;
    CombiData       data;
    
    fMCD12Q1Name = "/media/MEDIA/DataNCKH/MCD12Q1.A2011001.h27v06.051.2012264182800.hdf";

    findDimension(fMCD12Q1Name, &dimDataset, "Land_Cover_Type_1");
    dy =    (int) dimDataset[0];
    dx =    (int) dimDataset[1];
    
    // MCD12Q1
    land_cover_type1              = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type2              = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type3              = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type4              = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type5              = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type1_assessment   = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type2_assessment   = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type3_assessment   = (uint8 *)malloc(sizeof(uint8) * dx * dy);
    land_cover_type4_assessment   = (uint8 *)malloc(sizeof(uint8) * dx * dy);    
    land_cover_type5_assessment   = (uint8 *) malloc(sizeof(uint8) * dx * dy);
    land_cover_type_qc               = (uint8 *) malloc(sizeof(uint8) * dx * dy);
    land_cover_type1_secondary          = (uint8 *) malloc(sizeof(uint8) * dx * dy);
    land_cover_type1_secondary_percent        = (uint8 *) malloc(sizeof(uint8) * dx * dy);
    lc_property_1   = (uint8 *) malloc(sizeof(uint8) * dx * dy);
    lc_property_2   = (uint8 *) malloc(sizeof(uint8) * dx * dy);
    lc_property_3   = (uint8 *) malloc(sizeof(uint8) * dx * dy);    
    
    if (land_cover_type1 == NULL || land_cover_type2 == NULL || land_cover_type3 == NULL || land_cover_type4 == NULL || land_cover_type5 == NULL || land_cover_type1_assessment == NULL || land_cover_type2_assessment == NULL || land_cover_type3_assessment == NULL || land_cover_type4_assessment == NULL || land_cover_type5_assessment == NULL || land_cover_type_qc == NULL || land_cover_type1_secondary == NULL || land_cover_type1_secondary_percent == NULL || lc_property_1 == NULL || lc_property_2 == NULL || lc_property_3 == NULL ){
        printf ("Malloc Error");
        return (1);
    }
    
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_1", 0, land_cover_type1);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_2", 0, land_cover_type2);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_3", 0, land_cover_type3);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_4", 0, land_cover_type4);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_5", 0, land_cover_type5);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_1_Assessment", 0, land_cover_type1_assessment);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_2_Assessment", 0, land_cover_type2_assessment);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_3_Assessment", 0, land_cover_type3_assessment);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_4_Assessment", 0, land_cover_type4_assessment);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_5_Assessment", 0, land_cover_type5_assessment);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_QC", 0, land_cover_type_qc);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_1_Secondary", 0, land_cover_type1_secondary);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "Land_Cover_Type_1_Secondary_Percent", 0, land_cover_type1_secondary_percent);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "LC_Property_1", 0, lc_property_1);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "LC_Property_2", 0, lc_property_2);
    resHdf= getHDFMODPixelValueT (fMCD12Q1Name, "LC_Property_3", 0, lc_property_3);
        

    //string query;
    //stringstream query_stream;    

    /* Extract combination data
    foreach pixels in range
        print value to make maps
        update points table (land_cover value)
    */


    int sample, line;
    stringstream query_stream;
    for ( sample = top ; sample <= bottom ; sample++ ){
        for ( line = left ; line <= right ; line++ ){
            pos     = sample*dx + line;                   
            data.land_cover_type1               = (uint8)((long int)land_cover_type1[pos]);
            data.land_cover_type2               = (uint8)((long int)land_cover_type2[pos]);
            data.land_cover_type3               = (uint8)((long int)land_cover_type3[pos]);
            data.land_cover_type4               = (uint8)((long int)land_cover_type4[pos]);
            data.land_cover_type5               = (uint8)((long int)land_cover_type5[pos]);
            data.land_cover_type1_assessment    = (uint8)((long int)land_cover_type1_assessment[pos]);
            data.land_cover_type2_assessment    = (uint8)((long int)land_cover_type2_assessment[pos]);
            data.land_cover_type3_assessment    = (uint8)((long int)land_cover_type3_assessment[pos]);
            data.land_cover_type4_assessment    = (uint8)((long int)land_cover_type4_assessment[pos]);            
            data.land_cover_type5_assessment    = (uint8)((long int)land_cover_type5_assessment[pos]);
            data.land_cover_type_qc             = (uint8)((long int)land_cover_type_qc[pos]);
            data.land_cover_type1_secondary     = (uint8)((long int)land_cover_type1_secondary[pos]);
            data.land_cover_type1_secondary_percent     = (uint8)((long int)land_cover_type1_secondary_percent[pos]);
            data.lc_property_1          = (uint8)((long int)lc_property_1[pos]);
            data.lc_property_2          = (uint8)((long int)lc_property_2[pos]);
            data.lc_property_3          = (uint8)((long int)lc_property_3[pos]);
            
            printf ("%2d ", (int) data.land_cover_type1);

            query_stream.str(""); //clear
            query_stream << "UPDATE public.points SET land_cover=" << (int)data.land_cover_type1 << " WHERE line=" << line << " AND sample=" << sample << ";";
            postgresql_query( query_stream.str() );

        }
        printf("\n");
    }
   

    free (land_cover_type1);  
    free (land_cover_type2); 
    free (land_cover_type3); 
    free (land_cover_type4);
    free (land_cover_type5); 
    free (land_cover_type1_assessment); 
    free (land_cover_type2_assessment); 
    free (land_cover_type3_assessment); 
    free (land_cover_type4_assessment);     
    free (land_cover_type5_assessment); 
    free (land_cover_type_qc); 
    free (land_cover_type1_secondary); 
    free (land_cover_type1_secondary_percent); 
    free (lc_property_1); 
    free (lc_property_2); 
    free (lc_property_3); 

    time_t t2 = time(0);
    cout << "Total time processing : " << t2-t1 << endl;
    return 0;
}





// ******************************************************
// Read HDF file

GByte getHDFMODPixelValueT(char *fMODISName, char *bName, int bNum, VOIDP outRaster){

    int32 SDidentifier; // Scientific Data identifier

    if (( SDidentifier = SDstart(fMODISName, DFACC_READ)) == FAIL ){
        printf("ERROR: Unable to open inputfile '%s'\n", fMODISName);
        exit(1);}

    int32 nDatasets;
    int32 nGlobalAttrs;
    SDfileinfo(SDidentifier, &nDatasets, &nGlobalAttrs);

    int out_index= -1;
    int i;

    int32 DataSetId; // Data Set identifier
    char DataSetName[255]; // Name of the data set
    int32 DataSetDimensionNumber; // Number of dimensions in the data set; the maximum value is MAX_VAR_DIMS (or 32)
    int32 DataSetDimensionSizes[3]; // Array containing the size of each dimension in the data set; Y e X
    int32 DataSetDataType; // Data type for the data stored in the data set
    int32 DataSetNumAttrs; // Number of attributes for the data set

    for(i=0; i < nDatasets; i++){
        if ((DataSetId = SDselect(SDidentifier, i)) == FAIL) {
            printf ("ERROR: Unables select dataset #%d\n", i);
            exit(2);}

        SDgetinfo(DataSetId, DataSetName, &DataSetDimensionNumber, DataSetDimensionSizes,
                  &DataSetDataType, &DataSetNumAttrs);

        if( !strcmp(DataSetName, bName)){
            out_index = DataSetId;
            //printf ("getHDFMODPixelValueT()--> DS: '%s'\n", DataSetName);
            i= nDatasets;}
        if (SDendaccess(DataSetId) == FAIL) {
            printf ("ERROR: Unable to close SDS\n");
            exit(3);}
    }   // end for

    if(out_index== -1){
        printf ("getHDFMODPixelValueT()--> ERROR: No Data Set named '%s'\n", bName);
        exit(4);}

    int32 edges[DataSetDimensionNumber]; // Array specifying the number of values to be read along each dimension
    int32 start[DataSetDimensionNumber]; // Array specifying from where start read values

    if (DataSetDimensionNumber == 3){
        start[0] = bNum;
        start[1] = 0;
        start[2] = 0;

        edges[0] = 1;
        edges[1] = DataSetDimensionSizes[1];    // 2030 in MOD021km
        edges[2] = DataSetDimensionSizes[2];    // 1354 in MOD021km
    }
    else {
        start[0] = 0;
        start[1] = 0;

        edges[0] = DataSetDimensionSizes[0];    // 2030 in MOD021km
        edges[1] = DataSetDimensionSizes[1];    // 1354 in MOD021km

    }

    intn SDcode;
    if((SDcode = SDreaddata(out_index, start, NULL, edges, (VOIDP)outRaster)) == FAIL){
        printf ("getHDFMODPixelValueT()--> ERROR: Loading error \n");
        exit(7);}

    return 0;
}

// ******************************************************
// Read Dataset Dimension

void findDimension(char *nameInputFile, int **arrayDimension, char *nameBand){

    int32 SDidentifier; // Scientific Data identifier
    if(( SDidentifier = SDstart(nameInputFile, DFACC_READ)) == FAIL ){
        printf("ERROR: Unable to open inputfile '%s'\n", nameInputFile);
        exit(1);}

    int32 nDatasets;
    int32 nGlobalAttrs;
    SDfileinfo(SDidentifier, &nDatasets, &nGlobalAttrs); // get information about the Input Dataset

    int32 DataSetId; // SubDataset identifier
    char DataSetName[255]; // Name of the data set
    int32 DataSetDimensionNumber; // Number of dimensions in the data set; the maximum value is MAX_VAR_DIMS (or 32)
    int32 DataSetDimensionSizes[3]; // Array containing the size of each dimension in the data set; Y e X
    int32 DataSetDataType; // Data type for the data stored in the data set
    int32 DataSetNumAttrs; // Number of attributes for the data set

    int i;
    int n;
    for(i=0; i< nDatasets; i++){ // access to all SubDatasets

        if((DataSetId = SDselect(SDidentifier, i)) == FAIL){ // if it's not possible open 'i' SubDataset
            printf ("ERROR: Unables select SubDataset #%d\n", i);
            exit(2);
        }
        SDgetinfo(DataSetId, DataSetName, &DataSetDimensionNumber, DataSetDimensionSizes, &DataSetDataType, &DataSetNumAttrs);
        if(0 == strcmp(DataSetName, nameBand)){
            *arrayDimension= (int *)malloc(sizeof(int)*DataSetDimensionNumber);
            for(n=0; n<DataSetDimensionNumber; n++)
                (*arrayDimension)[n]= DataSetDimensionSizes[n];
        }
    } // end Subdatasets for cycle

}

// ******************************************************
// Read Factor (scale and offset)

void getFactor(char * inputFileName, char * band, char *attribute, float **valueToSearch){

    int32 SDidentifier; // Scientific Data identifier

    if (( SDidentifier = SDstart(inputFileName, DFACC_READ)) == FAIL ){
        printf("ERROR: Unable to open inputfile '%s'\n", inputFileName);
        exit(1);}
    int32 nDatasets;
    int32 nGlobalAttrs;
    SDfileinfo(SDidentifier, &nDatasets, &nGlobalAttrs);

    int32 subDataSetId; // Data Set identifier
    char DataSetName[255]; // Name of the data set
    int32 DataSetDimensionNumber; // Number of dimensions in the data set; the maximum value is MAX_VAR_DIMS (or 32)
    int32 DataSetDimensionSizes[3]; // Array containing the size of each dimension in the data set; Y e X
    int32 DataSetDataType; // Data type for the data stored in the data set
    int32 DataSetNumAttrs; // Number of attributes for the data set

    int32 out_index= -1;
    int i;
    for(i=0; i < nDatasets; i++){
        if ((subDataSetId = SDselect(SDidentifier, i)) == FAIL) {
            printf ("ERROR: Unables select dataset #%d\n", i);
            exit(2);}

        SDgetinfo(subDataSetId, DataSetName, &DataSetDimensionNumber, DataSetDimensionSizes,
                  &DataSetDataType, &DataSetNumAttrs);

        if( !strcmp(DataSetName, band)){
            out_index = subDataSetId;
            // printf ("getReflectanceFactor()--> DS: '%s'\n", DataSetName);
            i= nDatasets;}
        if (SDendaccess(subDataSetId) == FAIL) { // to end a 'SDcreate' or a 'SDselect'
            printf ("ERROR: Unable to close SDS\n");
            exit(3);}
    } // end for

    if(out_index== -1){
        printf("getReflectanceFactor()--> ERROR: No Data Set named '%s'\n", band);
        exit(4);}

    int32 attr_index, num_type, count;
    char attr_name[64];
    int k;

    attr_index = SDfindattr(out_index, attribute);
    if (attr_index == FAIL){
        printf ("ERROR: SDfindattr failed %s\n", attribute);
        exit(5);}

    // printf("getReflectanceFactor()--> Found '%s' attribute\n", attribute);

    // Get information about the file attribute
    if (SDattrinfo(out_index, attr_index, attr_name, &num_type, &count) == FAIL){
        printf ("ERROR: SDattrinfo failed\n");
        exit(6);
    }

    *valueToSearch= (float *)malloc(sizeof(float) * count);
    if(valueToSearch== NULL){
        printf("ERROR: can not malloc()\n");
        exit(7);
    }

    float buffer[count];
    // --------------------------------

    // Read the attribute data
    if(SDreadattr(out_index, attr_index, buffer) == FAIL){
        printf ("ERROR: SDreadattr failed\n");
        exit(8);}
    
    for(k=0; k<count; k++){
        (*valueToSearch)[k]= (buffer[k]);
    }

    // Terminate access to the array
    int status = SDendaccess(out_index);
    if (status == FAIL){
        printf ("ERROR: SDendaccess failed\n");
        exit(9);}

    // Terminate access to the SD interface and close the file
    status = SDend(SDidentifier); // to end a 'SDstart'
    if (status == FAIL){
        printf ("ERROR: SDend failed\n");
        exit(10);}
}
// ******************************************************


bool postgresql_query(string query)
{
    try
    {
        connection C("dbname=mymodisdb user=demo1usr password=demo1pwd hostaddr=127.0.0.1 port=5432");
        if ( !C.is_open() )
            return false;

        work W(C);
        W.exec( query );
        W.commit();

        C.disconnect ();
    }catch (const std::exception &e){
        cerr << e.what() << std::endl;
        return false;      
    }
    return true;
}